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Abstract 


A two-dimensional isothermal mechanistic model of an anode-supported solid oxide fuel cell was developed based on button-cell geometry. The 
model coupled the intricate interdependency among the ionic conduction, electronic conduction, gas transport, and the electrochemical reaction 
processes. All forms of polarizations were included. The molecular diffusion, Knudsen diffusion, as well as the simplified competitive adsorption 
and surface diffusion were also considered. An electric analogue circuit was used to determine the effective hydrogen diffusivity. The model results 
showed good agreement with the published experimental data in different H,—H,O mixtures without any other calibrations after the parameter 
estimation according to the experimental data in baseline operating condition. The distributions of species concentration and current density were 
predicted and the effects of cathode area, gas components, and anode thickness on the cell performance were studied. 


© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


Solid oxide fuel cells (SOFCs) are energy conversion devices, 
which produce electricity and heat directly from electrochemi- 
cal oxidation of the fuel. A PEN (positive electrolyte negative) 
structure is the main component of SOFC, which usually consists 
of two ceramic electrodes, anode layer and cathode layer, sepa- 
rated by a dense electrolyte layer. The electrodes are porous to 
facilitate the transport of fuel and oxidant from gas channel to the 
three phase boundaries where electrochemical reactions occur. 
The electrolyte is dense to keep the gases separated to form 
the oxygen concentration difference between the anode and the 
cathode. Oxygen ions are produced at the three phase bound- 
aries near the cathode/electrode interface, and are transported 
by a solid-state migration mechanism through the electrolyte to 
the anode/electrolyte interface, where oxygen ions react with 
fuel. Products molecules then transport back to the fuel channel 
through pores. 

The relative high operating temperature of SOFC can result 
in highly efficient conversion to power and heat for cogen- 
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eration or for a bottom cycle [1]. Some demonstrations have 
indicated that the SOFC hybrid systems can have a higher 
efficiency, and have both less pollutant emissions and lower 
greenhouse gas emissions than other alternatives [2]. With rising 
fuel prices and stricter emission control regulations, these capa- 
bilities make SOFCs even more attractive. However, SOFCs also 
possess some serious disadvantages in fabrication, sealing and 
operation, etc. A number of researchers are focusing on inter- 
mediate temperature solid oxide fuel cells (IT-SOFCs) between 
823 and 1073 K, allowing for a wider range of materials and 
more cost-effective SOFC fabrication [3,4]. Anode-supported 
SOFC is one of the attempts to IT-SOFCs, in which the anode 
is the support structure and is the thickest component in the 
PEN. The main advantage of the anode-supported SOFCs is 
the substantially lower ohmic resistance of the electrolyte, 
and consequently the lower operating temperature. However, 
it is generally accepted that the area specific resistance of IT- 
SOFCs may be larger than high-temperature SOFCs. This is due 
to activation, and perhaps concentration overpotential, which 
can often outweigh the ohmic contribution [5]. Thus, many 
efforts in SOFC technology development have been devoted 
to improve the electrode structure and electrolyte, aiming at 
obtaining a higher and more stable electrochemical performance 
[6,7]. 
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Nomenclature 

c concentration (mol m~?) 

D diffusion coefficient (m? s7!) 

F Faraday constant (96,485 C mol” !) 

AG electrochemical reaction activation energy 
(Jmol!) 

i current density (A m~?) 

ig exchange current density (A m~?) 

I total current density 

L thickness 

Mi molecular weight of species i 

m,n coefficient in ig calculation 

n mole fraction 

N mole flux (mol m~? s7!) 

p pressure (Pa) 

Q source term of charge balance equations (A m~*) 

T average pore size (wm) 

R gas constant (8.314 J mol! K7!) 

Ri source term of mass balance equations 
(kg m~? s7!) 

STPB TPB active area per unit volume (m? m7?) 

T temperature (K) 

V voltage (V) 

Xi mole fraction of species i 

Wi mass fraction of species i 


Greek letters 


a transfer coefficient 

E porosity 

n overpotential (V) 

O relative coverage of hydrogen on the TPB surface 
p density (kg m7?) 

o conductivity (S m7!) 
T tortuosity 

Subscripts 

an anode 

avg average 

ca cathode 

elec electronic 

ion ionic 


inter interface 


Subscripts 
bulk bulk phase 
eff effective 


Kn Knudsen 

MS Maxwell-Stefan 

surf surface phase 

TPB three phase boundary 


Since experimental studies on SOFC are expensive, time- 
consuming and labor-intensive, mathematical model is essential 
for SOFC study. Detailed mathematic model can help 
researchers to understand the complex, coupled reactions and 


transport processes within the cell, and can be used to study 
the effects of operating condition and electrode design parame- 
ters on the cell performance. It is of great value for mathematic 
model in reducing the time and cost in SOFC development if the 
model can predict the detailed processes of reaction and trans- 
port incorporated with the experimental data. Many researchers 
have studied these complex transport and reaction processes 
through either experiment or numerical modeling. However, due 
to a lack in understanding the intricacy transport processes, and 
the calculation complexity, it is still difficult to build a detailed 
mechanistic model. Traditional models usually need to adjust 
many parameters to fit the experimental data, and when the 
operating conditions changed, the model parameters have to 
be adjusted again to fit the experimental data. Especially for 
anode-supported SOFC, the existing models in open literature 
usually calculate the concentration polarization by using a semi- 
empirical formulation [3,4]. The bulk diffusion considering the 
molecular and Knudsen diffusion are treated as the main sources 
of diffusion resistances. In this case, the tortuosities have to 
be tuned in the range of 10-17 to generate relatively higher 
anodic concentration overpotentials, then to fit to the experimen- 
tal polarization curves. However, the modern anode tortuosities 
measured in experiments are mainly in the range of 2.5-3.5 [8]. 
This suggests the tuned tortuosity values are not consistent with 
the fact that the tortuosity is only a physical property param- 
eter. Williford et al. [8] developed a semi-empirical model of 
anodic concentration polarization considering the competitive 
adsorption and surface diffusion. Although the model applica- 
bility should be further improved and the transport mechanisms 
should be further clarified through more detailed experiments, 
it still suggested that the diffusion resistance of modern porous 
ceramic anode materials is not only in the bulk of the material, 
but also is near the electrode/electrolyte interface. Models incor- 
porating such concepts can predict the limiting current behavior 
and the performance envelopes more accurately. 

In this paper, a two-dimensional isothermal mechanistic 
model of an anode-supported solid oxide fuel cell to describe 
the intricate interdependency among the ionic conduction, elec- 
tronic conduction, gas transport and electrochemical processes is 
presented. The molecular diffusion, Knudsen diffusion, compet- 
itive adsorption near the reactive sites, and the surface diffusion 
are considered to escape from the unreasonable adjustment of 
tortuosities in various H2-H20 mixtures. To improve the model 
generality, the general form of Butler—-Volmer equation was used 
to solve the cell current densities and all forms of polarizations 
numerically. The exchange current densities were formulated as 
a general formulation. According to the polarization curve in a 
predefined base case, some parameters, which are difficult to be 
achieved from experiments or from published literatures were 
estimated. 


2. Model development 
2.1. Model assumptions 


Consider a typical anode-supported solid oxide fuel cell, as 
shown in Fig. 1. 
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Electrolyte layer 
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O electronic conductor 
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Fig. 1. Schematic diagram of anode-supported SOFC button cell. 


Model assumptions are listed as follows: 


(1) Steady state conditions assumption is applied in this model; 

(2) The reactant gas mixtures are approximated as ideal gases, 
thus gas mixture physical properties such as specific vis- 
cosity, density, etc., can be easily estimated according to the 
mixed gas composition; 

(3) Since the button cell is relative small, the temperature vari- 
ation across the cell is small, thus model is assumed as 
isothermal and the physical properties are evaluated at the 
average cell temperature; 

(4) The reaction active sites are uniformly distributed in the 
electrode. The two conducting phases, electronic conduct- 
ing phase and ionic conducting phase, are considered to be 
continuous and homogeneous. Since the material properties 
of both electronic conductor and ionic conductor are only 
represented by model parameters, the model itself could be 
assumed to be not material-specific; 

(5) The pressure gradient in the porous flow is neglected. 


Furthermore, the boundary conditions of concentration 
and potential is assumed uniform at the electrode/gas channel 
interface. 


Z4 


OQ 
OQ 
OQ. 


OCA eee chamber 


ca/air chamber 


ca/elec 


an/elec 


2.2. Model geometry 


The model geometry is based on a button cell. Due to the 
cell symmetry conditions, the two-dimensional axial symmetry 
coordinate is adopted as that illustrated in Fig. 2. 


2.3. Charge balance 


2.3.1. General formulation of charge balance governing 
equation 

The charge balance control domains include the anode, the 
cathode and the electrolyte. At the electrodes, electrons and oxy- 
gen ions are served as conductive particles, and at the electrolyte, 
the conductive particles are only oxygen ions. The charge bal- 
ance can be formulated by the point form of Ohm’s law in a 
stationary coordinate system: 


-V -(6VV)= Qj (1) 


The current source term Q; can be calculated by the equation 
below. 


Q j = +itrans STPB (2) 


elec/air chamber 


Da 


aQ So 


symmetry 


insulation 


r 


Fig. 2. Schematic diagram of model geometry. 
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The choice of the sign before the term on the right side of 
equation depends on the conductor phase type and calcula- 
tion domain. For ionic charge balance at the anode, the current 
source term should be positive as the anode is served as the ion 
sources and the ions carries negative electric charge, while, for 
electronic charge balance at the anode, the current source term 
should be negative to keep the charge balance. Similar to the 
anode, the cathode is served as ion sinks and electron sources, 
thus, the current source term is negative when considering the 
ionic charge balance and is positive when considering electronic 
charge balance. 

Stpp in above equation is the electrochemical reaction active 
area per unit volume, this parameter has been formulated by 
some researchers by using the particle coordination number, 
together with percolation theory as follows [9-11]. However, to 
simplify the calculation, this parameter was kept as a constant 
in this paper. 

The itrans is the local charge transfer current density and can 
be calculated by generalized Butler-Volmer equation: 


i i (= exp (« ne Fn 
trans = 10 T 
Creact RT 


ae exp ( a ane) (3) 


Cprod RT 


where ig is the exchange current density; a the charge transfer 
coefficient; F the Faraday constant; R the gas constant; Creact; 
ae Cprod, Gori are the reactant and product concentrations at 
the reaction active sites and the electrode/channel interface; 7 
the local overpotential, defined as follows [11]: 


n = Velec — Vion — Viet (4) 


where Vref is the relative potential difference between the 
electronic and ionic conductors at the reference state. In this 
paper, the open-circuit state was chosen as the reference 
State. 


2.3.2. Charge balance at the cathode layer 

The cathode usually gives large contribution to potential loss. 
For conventional cathode material, electrochemical reaction 
transfers the current from ionic conductor to electronic con- 
ductor. The voltage loss mainly consists of Ohmic overpotential 
in the electronic and ionic conductors, activation overpotential 
for electron transfer reaction, and concentration overpotential, 
due to mass transport resistance. To deal with the concentra- 
tion overpotential, the gas species concentration differences 
between the reactive sites and electrode/air channel interface are 
considered. 


girs 
Creact,ca _ ©Qp,ca TPB _ „bulk 5 
I — bulk ” Cprod,ca = Cprod,ca (5) 
Creact,ca COy,ca 
where Coe and ok are the oxygen concentrations at the 


reaction active site and at the electrode/gas chamber interface. 
The exchange current density is defined as the current density 

of the charge-transfer reaction at the dynamic equilibrium poten- 

tial. Under this circumstance, the forward and reverse current 


densities are equal to ip. Thus ig is a measure of the electro- 
catalytic activity at the electrode/electrolyte interface or the TPB 
for a given electrochemical reaction [12]. The cathodic exchange 
current density can be formulated as [13]: 


ye 


(Po2,ca/ PC 
A 2,ca O2,ca ip (6) 


10,ca = lo, ä 
1 + (P02,ca/ Pb,,ca) 


where il ca 18 used as an empirical parameter to fit experimental 
observation; po, the local oxygen partial pressure; POs, ca IS 
formulated as the Arrhenius form as below: 


=) (7) 


Ponca = Ao, exp ( RT 


where Ao, = 4.9F8 bar, and Eo, = 200 kJ mol~!, 
By setting Vref,an to zero, the cathode reference potential 
Vref,ca Will be the open circuit voltage: 


0.5 
RT H2, 
Vief,ca = E? In Tones 
2F PH)0,a 


(8) 


where Z? is the ideal Nernst potential and can be calculated as: 


AG. 
nNeF 


E? = 


(9) 


where AG° is the change of standard-state Gibbs free energy 
between the products and the reactants. 

However, calculated open circuit voltages usually are a little 
higher than the experimental values. In this paper, this difference 
was treated as the so called “leak overpotential” denoted by neak- 

Then, the ionic charge balance at the cathode can be written 
as: 


=V. (Gion,caV Vion, ca) 


TPB 
i S. Co, ane F(Velec,ca aa Vion,ca = Vief,ca) 
= —10,caSTPB,ca 
cbulk RT 


Oz 
ls eF Valec ca T Vion ca Vie ca 
xex ( Q = a@)ne F (Va s ; f, ’)) (10) 


Similar to the ionic charge balance, the electronic charge balance 
at the cathode can be written as: 
-=V .- (Gelec,ca V Valec,ca) 


TPB 
i S Co, ex ane F (Velec,ca = Vion,ca = Vref,ca) 
= l0,ca9TPB,ca Take 
cbulk RT 


(07) 
1— eF Velec,ca — Vion,ca — Vre ca 
X exp ( (1 — a)ne F( ma , f, )) (11) 


2.3.3. Charge balance at the anode layer 

Similar to cathode, for conventional anode material, Ni/YSZ 
mixture, electrochemical reaction transfers the current from 
ionic conductor YSZ to electronic conductor Ni. 
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According to the general governing equation of charge bal- 
ance, the ionic current balance at the anode is given as: 


=V. (Gion,an V Vion,an) 


TPB 
i S (= exp ane F(Vetec,an K Vion,an a Vref,an) 
= —l0,an9TPB,an bulk 
cH RT 


co ad = )Ne F(Velec,an > Vion,an B Vief,an) 
bulk ex RT 
H0 
(12) 
The calculation of hydrogen species at the TPB active sites will 
be discussed in Section 2.4. 
The exchange current density of hydrogen electrochemical 
oxidation was proposed in this paper. 


i : CHD — 120000 
10,an = id an (z P ) exp (=) (Posen)? (13) 
ase, H2 


where i en is an adjustable parameter to fitting the experimental 
data in base case. Cpase,H, the hydrogen concentration in the base 
case. 

The po,,an is the water vapor partial pressure at the anode 
due to hydrogen oxidation equilibrium and could be calculated 
as: 


2 
H20 
Pozan = (=) (14) 
Keg,H> PH> 


where Keg,H, is the equilibrium constant of hydrogen oxidation 
reaction. 

Similar to the ionic charge balance, the electronic charge 
balance at the anode can be written as: 


—V + (Oelec,an V Velec,an) = —itrans,an STPB,an (15) 


2.3.4. Charge balance at the electrolyte layer 

For the electrolyte layer, there is only ion conductor. And 
there are no current sources or sinks. Consequently, only Ohmic 
polarization exists in the electrolyte and the charge balance gov- 
erning equation can be written as: 


ff 
—V. (ion electrolyte V Vion,electrolyte) =0 (16) 
where O elecnoiyte is the effective conductivity of the YSZ 


electrolyte. 
2.4. Mass balance 


2.4.1. General formulation of mass balance governing 
equation 

Diffusion inside the porous media is important for fuel cell 
reactions. Within a pore, three fundamentally different types of 
diffusion mechanisms are distinguished [14]. 


(1) Free molecular diffusion is significant for large pore size and 
high pressures. In this case, molecule-molecule collisions 
dominate over molecule-wall collisions. 

(2) Knudsen diffusion becomes predominant when the mean- 
free path of the molecular species is much larger than the 


pore diameter and hence molecule-wall collisions become 
important. 

(3) Surface diffusion of adsorbed molecular species along the 
pore wall surface becomes dominant for micro pores and 
for strongly adsorbed species. At high currents and high 
fuel utilization, the competitive adsorption and surface dif- 
fusion of hydrogen may be responsible to the high interface 
diffusion resistance [8]. 


The steady state mass balance can be formulated: 
V- ji= Ri (17) 


where R is the source term of mass balance, and can be formu- 
lated as the reaction rate of the electrochemical and chemical 
reactions; j the molecular mass flux and can be formulated as 
a function of a diffusion driving force. For species i, accord- 
ing to the generalized formulation of Fick’s law, the governing 
equations of mass balance could be expressed as: 


N 
-V . (nS 0H) = R; (18) 
k=1 


where w; is the mass fraction of species i; dg the generalized driv- 
ing force. Considering the driving forces through concentration 
and pressure, the d can be formulated as: 


dy = V xk (19) 


pof is the effective multi-component diffusivity between com- 
ponents i and k. In order to solve Eq. (18), the multi-component 
diffusivity must be specified. 

Free molecular and Knudsen diffusion mechanisms occur 
together and it is prudent to take both of them into account 
rather than assume that one or the other mechanism is con- 
trolling. It is generally agreed that the dusty gas model is the 
most convenient approach to describe the combined bulk and 
Knudsen diffusion [15]. The principle of the dusty gas model is 
based on the Maxwell—Stefan diffusion by considering the pore 
wall as consisting of giant molecules uniformly distributed in 
space. In order to account for the tortuous path of the molecule 
rather than along the radial direction and the porosity of the elec- 
trode for the fact that diffusion occurs only in the pore space, 
an effective Maxwell—Stefan diffusion coefficient and Knudsen 
diffusion coefficient can be formulated as: 


E 

pie = — Dike (20) 
E 

pret = = pe (21) 
T 


where, £ denotes the electrode porosity and t denotes the elec- 
trode tortuosity. 

Then, the bulk effective multi-component diffusivities can be 
calculated by: 


1 

bulk, eff 

Dv = (22) 
ik 1/DMSef ip 

For species O2, N2 and H20, the effective diffusivities are just 

the bulk effective multi-component diffusivities. For hydrogen 
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Molecular free diffusion+ 
Knudsen diffusion 


Total diffusion flux 


Bulk (electrode/fuel 
chamber interface) 
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electrode) 


1+A1(Al<<1) 


Surface diffusion 


Molecular free diffusion+ 
Knudsen diffusion 


Fig. 3. Electric analogue circuit of H2 diffusion. 


species, to achieve the H2 concentration at the TPB, cH the 
competitive adsorption and surface diffusion T Da were 
taken in to account. 

The electric analogue circuit of hydrogen diffusion is 
assumed as shown in Fig. 3. 

As the figure shows, the diffusion flux from fuel chan- 
nel/electrode interface to a certain place within the porous 
electrode is equal to the diffusion flux from this certain place 
to the three phase boundary. Here, / is the distance from 
anode/electrolyte interface to a certain position within the elec- 
trode and AŻ the distance from a certain position within the 
electrode to the three phase boundaries. A/is much shorter than /, 
and the effective diffusivity D can be formulated as following 
Eq. [8]: 


= (Dy eff)? om? (23) 


where @ is the relative coverage of hydrogen on the TPB surface, 
pet the surface diffusion coefficient and can be formulated as 


[8]: 


(Det) KON 


1-0 


pet = = 


(24) 


where Diy is the surface diffusion coefficient of hydrogen as 


zero coverage is approached, and Di is that as full coverage is 
approached at very low relative coverage and high temperatures, 
DE") may approach the bulk gas value due to the high avail- 
ability of vacant surface sites for diffusion by atomic hopping 
mechanisms. Thus, to simplify the calculation, it is reasonable 
to assume the DEG is equal to Dy ett Also, we assumed 
that the relative coverage of hydrogen is approximately equal to 
mole fraction of hydrogen in this paper. The ratio of D“ i 


and Dat remains to be an adjustable parameter in the model. 


2.4.2. Mass balance at the cathode 

According to the Faraday’s law, the relationship between the 
mass balance source term and the current source term at the 
cathode can be built: 


itrans,O. STPBMo, 


Rop 4F 


(25) 


The reaction rate of nitrogen at the cathode is zero and the mass 
fraction of nitrogen can be calculated by: 


WN, =1- WO, (26) 


2.4.3. Mass balance at the anode 

At the anode, the reactions include the electrochemical 
reactions of hydrogen and carbon monoxide and the shifting 
reactions. 

Similar to the cathode, the relationship between the mass 
balance source term and the current source term for the electro- 
chemical reactions at the anode can be built: 


itrans,an STPB Mp, 


2F 27) 


Ry, = 


—itrans,an STPB MH20 
2F 


Rmo = (28) 
From the mass balance at the anode, the species concentration 
at the electrode, c;, can be solved. 


2.5. Boundary conditions 


In order to solve the coupled partial differential equation sys- 
tems of charge and mass balance, the boundary conditions of all 
outer interfaces are specified as in Table 1. 

In Table 1 the difference between Veell.ca and Veell,an Should 
be the cell voltage used in the calculation. Here, we choose 
Veellan = 0 V. Wi, buik in the table is the mass fraction of species 
i in the fuel/air channel. The boundary condition “insulation” 
denotes no flux through the boundary. the boundary condition 
“symmetry” denotes that the partial derivate of the parameter at 
the boundary is zero. 


2.6. Solution algorithm 


By given a specified value of cell voltage, the parameters 
(e.g. current density, species concentration) distributions can be 
achieved. The calculation was done by using the finite element 
commercial software COMSOL MULTIPHSICS®. 
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Table 1 
Boundary conditions 


Ionic charge balance 


Electronic charge balance Mass balance 


O22 calair chamber Insulation Veell,ca WO , bulk; WN, bulk 
O82 decisis chamber Insulation Insulation Insulation 

OQ ealelec Continuity Insulation Insulation 

OMan/elec Continuity Insulation Insulation 

OL2an/fuel chamber Insulation Veell,an WH), bulk; WH>0,bulk 
IN symmetry Symmetry Symmetry Symmetry 

3NR others Insulation Insulation Insulation 

Table 2 

PEN geometries and properties parameters 

Layer Cathode Electrolyte Anode 

Radius, R (cm) 0.8 1.3 1.3 

Thickness, d (um) 60 10 1100 

Porosity, € 0.35 - 0.35 

Tortuosity, T 3.5 - 35 

Average pore radius, rp (Įm) 0.5 - 0.5 

Ionic conductor conductivity, Cion (S m7!) 3.34 x 104 exp(—10,300/T) 3.34 x 104 exp(—10,300/T) 3.34 x 10* exp(—10,300/7) 
Electronic conductor conductivity, Celec (S m7!) seit exp (=) - 2x 10° 


According to the current density distribution, the average 
current density can be calculated by: 


Trote 1 fk 
Iavg = a = I 2rilocal dr 
0 


(29) 


where Rinter is the radius of the interface we concerned; ioca] the 
local current density. 

In order to generate the full polarization curve, the calculation 
should be done over a range of cell voltages to calculate the 
corresponding average current density. 


3. Results and discussions 
3.1. Parameter settings and model validations 


Some of the model input parameters are listed in Table 2. 

The operation conditions in the base case are listed in Table 3. 
Since some of the parameters could not be achieved directly from 
experiments or published literature, the tuning processes have to 
be involved in the study. One of the adjustable parameter, leak 
overpotential, could be estimated as the difference between the 
experimental and calculated open circuit voltage using Eq. (9). 
This parameter was kept as a constant 0.05 for all of the calcula- 
tions. The adjustable parameters ib ane i i ca and Srpg were varied 
in order to bring the model results into closer agreement with 


Table 3 

Operating conditions in the base case 

Parameters Value 
Pressure, p (Pa) 101,325 


Temperature, T (°C) 800 
Fuel composition 85% Hz and 15% H20 
Oxidant composition 21% Oz and 79% No. 


the experimental data in the base case. It should be noted that 
the adjustable parameter De / pat , was varied in order to 
assure the agreement of the modeling results and experimen- 
tal data in the case when the fuel composition is 20% H2 and 
80% H20. Then the model could be validated further in some 
other cases without changing these parameters. The values of 
adjustable parameters are summarized in Table 4. 

Then, by changing the fuel composition and cell voltage, the 
cell performances obtained from the model for various compo- 
sitions are compared with the experimental data at 800°C with 
air as the oxidant as shown in Fig. 4: 

The results showed that the calculated polarization curves 
agreed very well with the experimental data in most cases except 
that in the case of pure hydrogen and the case of 50% H2, 50% 
H20. For the 50% H2—-50% H20 fuel mixtures, the model results 
overestimate the cell performance when the current density is 
relatively high. And for the pure hydrogen, the model results 
slightly underestimate cell performance in the overall extent 
of current density. In fact, the simplification to the competi- 
tive adsorption and surface diffusion procedure might be the 
important source of the derivation. Owing to the consideration 
of detailed processes of adsorption and surface diffusion, the 
model can be used to predict the cell performances in a large 
extent of current density while also in a large fuel composition 


Table 4 

Adjustable model parameters 

Parameter Value 
Leak overpotential, leak (V) 0.05 
Adjustable parameter, ih, an (A m~?) 1.079E4 
Adjustable parameter, ib, ca (A m~?) 2E3 
TPB area per unit volume, Srpg (m? m~?) 3.225E5 


Adjustable parameter, Da / oe i 8 
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Fig. 4. Results of parameter estimation and model validation. 


extent without any other calibrations after the parameter estima- 
tion according to the experimental J-V curve in the base case. 
The results suggested that with the amount of H2 increasing 
in the fuel, the open circuit voltage as well as the maximum 
current density decreased. When H2 concentration is low, sub- 
stantial concentration polarization was present as evidenced by 
the observation of a limiting current density. 

For making detailed comparison, Fig. 5 shows the modeling 
results without the modifications to the hydrogen concentration 
at the three phase boundaries while using the same parametric set 
with baseline condition except that the tortuosity is treated as an 
adjustable parameter to assure the agreement of the modeling 
results and experimental data. We can see that when the fuel 
composition is 85% H2 and 15% H20, the value of tortuosity is 
3.5 and this value agreed well with the published experimental 
data. However, when the fuel composition is 20% H3 and 80% 
H20, the model could not predict the experimental data well 
until the tortuosity in the model was adjusted to 14. 


3.2. Species concentration distribution in the base case 


The distributions of the H2, H2O concentrations at the anode 
and the O2 concentration at the cathode in the base case are 
shown in Fig. 6. 
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Fig. 5. Modeling results by tuning tortuosity. 
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Fig. 6. Surface plots of H2, H20, O2 concentration. 


The H, concentration reduces fast along the z direction in the 
region opposite to the cathode and reduces more gradually in 
the region outside of cathode radius, and this un-uniformity is 
due to the asymmetric electrode of the button-cell. Contradict to 
Hp distribution, the H2O concentration rises along Z direction. 
The O2 concentration distribution remains relatively uniform 
throughout the cathode. This suggests that the O2 is sufficient 
even just using air as oxidant for electrochemical reaction and 
that the cathodic concentration overpotential can be neglected 
compared to relatively large anodic concentration overpotential. 


3.3. Current densities distribution 


Fig. 7 gives the electronic current density and ionic cur- 
rent density distributions at the anode, cathode and electrolyte 
along the cell thickness direction (r= 0). As shown in the figure, 
most of the electrons are generated near the electrode/electrolyte 
interface where ions are eliminated. The calculations sug- 
gest that the electrochemical reactions mainly occur at the 
electrode/electrolyte interface. It should be noted that with 
decreasing of H2 content in the fuel, the reaction occurs at a 
position further from the anode/electrolyte interface. In fact, 
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Fig. 7. Distributions of electronic and ionic current density along the button-cell 
thickness direction (r=0). 
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Fig. 8. Effects of cathode area on the cell performance for (a) fuel composition: 85% Hz and 15% H20; (b) fuel composition: 20% Hz and 80% H20. (anode area: 


5.31 cm?). 


the reactions always occur at the sites where the total resis- 
tances, including mass transport resistance, electrochemical 
activation resistance and also the charge transfer resistance, 
reach the minimum [16,17]. Lower H2 concentration will lead 
to larger diffusion resistance, so at this time, it is much more 
difficult for the reactant gas to approach the anode/electrolyte 
interface. 


3.4. Effects of cathode area on cell performance 


Fig. 8 shows cell polarization curves as well as power density 
versus current density traces, for cells with cathode areas of 0.5, 
2, 3.5 cm? at 800°C. Both the current density and power density 
are determined on the basis of the cathode area. 

As shown in Fig. 8(a), the fuel composition is 85% H2 and 
15% H20. There was no obvious limiting current density even 
at the highest current density. The maximum power density is 
approximately 16,500 W m~? with the 0.5 cm? cathode, while 
that for the cell with the 2cm? cathode is ~14,300 Wcm~?. 
The results show that the difference is 2200 Wm, with the 
relative difference is ~13.3%. When the cathode area increase 
from 2 to 3.5 cm?, the difference of maximum power density is 
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about only 100 W cm~*. This difference can almost be neglected 
even the cathode area is doubled. Similar to this case, as shown 
in Fig. 8(b), the fuel composition is 20% H2, 80% H20.The 
maximum power density is 6770 W m~? with 0.5 cm~?, while 
that for the cell with the 2 cm? cathode is about ~4930 W m7?. 
The difference is about 1840 Wcm~?, and the relative differ- 
ence is almost ~27.2%. When the cathode area increase from 2 
to 3.5 cm?, the difference is about 180 W m~, and the relative 
difference is about ~3.6%. Thus, a preliminary conclusion is 
that the cell performance will be better with decreasing cath- 
ode area while keeping anode area as a fixed valued. The effect 
degree of cathode area will be higher with the smaller cathode 
or with the lower H? content fuel mixtures. This conclusion is 
qualitatively constant with the analysis in the published liter- 
ature [4]. The source of potential differences in the reported 
power densities due to asymmetric electrodes lies principally 
in the differences in the anodic concentration polarizations. It 
suggests that the relatively larger anode/fuel gas interface with 
smaller active anode/electrolyte interface could lead to the lower 
gas transport resistance. In order to predict the anodic concen- 
tration polarization correctly, the cathode area should not be too 
small in the real experiment. 
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Fig. 9. Effects of anode thickness on the cell performance for (a) fuel composition: 85% Hz and 15% H20; (b) fuel composition: 20% Hz and 80% H20. 
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3.5. Effects of anode thickness with different fuel 
composition 


Fig. 9 presents the effects of anode thickness on the cell 
performance. As shown in Fig. 9(a), when the fuel composi- 
tion is 85% H and 15% H2O mixtures, the cell performance 
is slightly affected by the anode thickness and there was no 
obvious limiting current density even when the relatively thick 
anode (1800 pm) was adopted. Since the effects of anode thick- 
ness are mainly on the anodic concentration overpotential, it 
suggests that the anodic concentration overpotential is not obvi- 
ous in this case. However, if the reactant gas concentration in 
the fuel is low, as Fig. 9(b) shows, the cell limiting current den- 
sity decrease fast from 38,400 to 11,200 A m7? when the anode 
thickness increase from 300 to 1800 um. Thus, in actual oper- 
ation of SOFC, when the reactant concentration is low or fuel 
utilization is high, especially in the fuel channel downstream 
area, gas diffusion through a thick and porous anode will be one 
of the performance-limiting factors. In this case, a thinner anode 
layer should be adopted to reduce the large anodic concentration 
overpotential. 


4. Conclusions 


A two-dimensional isothermal mathematical model of an 
anode-supported solid oxide fuel cell was developed based 
on a button-cell geometry. The model coupled the intricate 
interdependency among the ionic conduction, electronic con- 
duction, gas transport, and electrochemical processes. The 
model has been applied to generate polarization curves and 
parameter distributions for different cell configurations and 
specified operating conditions. All forms of the polarizations 
were considered, including activation polarization, Ohmic polar- 
ization, concentration polarization, and leak overpotential. For 
the modifications to the mass transport calculation, the molec- 
ular diffusion, Knudsen diffusion, as well as the simplified 
competitive adsorption and surface diffusion processes were 
taken into account. The validation results showed that the model 
agrees well with the published experimental data in different 
H»—H20 fuel mixtures without any other calibrations after the 
parameter estimation according to the experimental data in base- 
line operating condition. 

The distributions of H2, H2O and Op at the electrode sug- 
gest that the O2 is sufficient even by using air as oxidant. The 
cathodic concentration overpotential can be neglected compared 
to relatively large concentration difference at the anode if we 
want to simplify the calculation. The results of electronic/ionic 


current densities distributions suggest that the electrochemical 
reactions mainly occur near the electrode/electrolyte interface 
and with decreasing H» content in the fuel, the reaction occurs 
at a position further from the anode/electrolyte interface. The 
cell performance will be better with a smaller cathode while 
keeping the anode area as a fixed value. The cathode area effect 
will be higher with the smaller cathode or with the lower H2 con- 
tent fuel mixtures. The effects of anode thickness is not obvious 
when the H3 content in the fuel mixtures is relatively high, while 
that is very significant if the H2 content in the fuel mixtures is 
low. This situation usually occurred at the downstream of fuel 
channel especially when the fuel utilization is high. 

Finally, it should be noted that it is necessary to improve 
the competitive adsorption and surface diffusion mechanism to 
predict the detailed mass transport processes within the porous 
anode. In addition, in future work, the model will be extended 
to study the cell performance with multi-component fuel, which 
may be more practical in actual SOFC applications. 
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